%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Autler-Townes effect eigenvalues & Bloch
%
% Stationary optical Bloch equations for a Lambda system,
%  n = [22 33 12 21 13 31 23 32],
%  ---	 |3>
%    \
%    --- |2>
%    /
%  ---	 |1>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch
    figure(1);          clf;
    col='kbgrmcy';

    Gam=1;
    omg=[0:.01:1]*Gam;

    % eigenvalues

    for kk=1:length(omg)
        Omg=omg(kk);
        Eigen = eig( [0 Omg/2 ; Omg/2 -1i*Gam/2] );
        [Shift,I] = sort(real(Eigen));
        Width = -sort(imag(Eigen(I)));
        E1(kk)=Shift(1);        W1(kk)=Width(1);
        E2(kk)=Shift(2);        W2(kk)=Width(2);
    end

    subplot(331),plot(omg/Gam,E1,omg/Gam,E2,'LineWidth',1)
    set(gca,'FontSize',16,'LineWidth',1)
    xlabel('\fontsize{16}\Omega/\Gamma')
    ylabel('\fontsize{16}\Delta{\itE}/\Gamma')
    text(-.3,.5,'\fontsize{16}(a)')
    subplot(332),plot(omg/Gam,W1,omg/Gam,W2,'LineWidth',1)
    set(gca,'FontSize',16,'LineWidth',1)
    xlabel('\fontsize{16}\Omega/\Gamma')
    ylabel('\fontsize{16}\Gamma_{eff}/\Gamma')
    text(-.3,.6,'\fontsize{16}(b)')

    
    % Bloch %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    G12=		       1;                                                   % decay of populations
    G23=		      .1;
    G13=		     .01;
    O12=   [0:.02:1.5]*G12;                                                 % Rabi frequencies
    o23=	      .1*G23;
    D23=[-1:.01:1]*7*G23;                                                   % detunings

    for kk=1:length(O12)
        o12=O12(kk);
        for ll=1:length(D23)
            d12=          0;
            d23=	D23(ll);
            n = Matrix_K(G12,G23,G13,o12,o23,d12,d23);
            nt(kk,ll)=real(n(2));
        end
    end

    subplot(333),mesh(D23/G12,O12/G12,nt)
    set(gca,'FontSize',16,'LineWidth',1.3)
    xlabel('\fontsize{16}\Delta_{23}/\Gamma_{23}');
    ylabel('\fontsize{16}\Omega_{12}/\Gamma_{12}');
    zlabel('\fontsize{16}{\it\rho_{33}}');
    text(-1.7,1,4e-4,'\fontsize{16}(c)')
    shading flat
    axis tight
    view([30 70])
end


function n = Matrix_K(G12,G23,G13,o12,o23,d12,d23)
    I=          1i/2;
    g12=         G12;                                                       % decay of coherences
    g13=     G13+G23;
    g23= G12+G13+G23;
    d13=     d12+d23;
    L12=g12/2+1i*d12;
    L13=g13/2+1i*d13;
    L23=g23/2+1i*d23;

    b=[0 0 -I*o12 I*o12 0 0 0 0]';
    
    M=[
       -G12     G23		-I*o12	I*o12	0		0		I*o23	-I*o23
        0		-G13-G23 0		0		0		0		-I*o23	I*o23
        -1i*o12	-I*o12	-L12	0		I*o23	0		0		0
        1i*o12	I*o12	0		-L12'	0		-I*o23	0		0
        0		0		I*o23	0		-L13	0		-I*o12	0
        0		0		0		-I*o23	0		-L13'	0		I*o12
        I*o23	-I*o23	0		0		-I*o12	0		-L23	0
        -I*o23	I*o23	0		0		0		I*o12	0		-L23'
    ];

    n = -M\b;
end